High-accuracy waveforms for binary black hole inspiral, merger, and ringdown 
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The first spectral numerical simulations of 16 orbits, merger, and ringdown of an equal-mass 
nonspinning binary black hole system are presented. Gravitational waveforms from these simulations 
have accumulated numerical phase errors through ringdown of < OT radian when measured from the 
beginning of the simulation, and < 0.02 radian when waveforms are time and phase shifted to agree 
at the peak amplitude. The waveform seen by an observer at infinity is determined from waveforms 
computed at finite radii by an extrapolation process accurate to < 0.01 radian in phase. The 
phase difference between this waveform at infinity and the waveform measured at a finite radius of 
r = lOOM is about half a radian. The ratio of final mass to initial mass is Mf/M = 0.95162±0.00002, 
and the final black hole spin is Sf/M^ = 0.68646 ± 0.00004. 
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I. INTRODUCTION 

Beginning with the groundbreaking binary black hole 
evolutions of Pretorius [1] and the development of the 
moving puncture method [2, 3], it has recently become 
possible to solve Einstein's equations numerically for the 
inspiral, merger, and ringdown of two black holes in a 
binary orbit. Already these simulations have provided 
tests of post-Newtonian approximations [4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14], have allowed initial exploration of 
the orbital dynamics of spinning binaries [15, 16, 17, 18, 
19, 20], have determined the recoil velocity of the final 
black hole when the masses are unequal [21, 22, 23, 24], 
and have led to the discovery of dramatically large recoil 
velocity from certain spin configurations [18, 25, 26, 27, 
28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. 

Waveforms from these numerical simulations are im- 
portant for gravitational-wave detectors such as LIGO 
and LISA. This is not only because detected waveforms 
can be compared with numerical models to measure as- 
trophysical properties of the sources of gravitational radi- 
ation, but also because the detection probability itself can 
be increased via the technique of matched filtering [38] , in 
which noisy data are convolved with numerical templates 
to enhance the signal. 

However, binary black hole simulations are time con- 
suming: a single simulation following approximately 10 
orbits, merger, and ringdown typically requires a few 
weeks of runtime on approximately 50 or 100 processors 
of a parallel supercomputer, and typically such a simula- 
tion produces waveforms of only modest accuracy. This 
large computational expense precludes, for example, pro- 
ducing a full template bank of numerical waveforms cov- 
ering the entire parameter space of black hole masses and 
spins. Hence there has been much interest in construction 
of phenomenological analytical waveforms [7, 39, 40, 41] 
that can be computed quickly and are calibrated by a 
small number of numerical simulations. While the accu- 
racy of typical simulations is sufficient for creating LIGO 



detection templates, it is most likely inadequate for LIGO 
parameter estimation and is far from what is required for 
LISA data analysis [42] . 

One approach to increasing the accuracy and effi- 
ciency of simulations is to adopt more efficient numer- 
ical methods. In particular, a class of numerical tech- 
niques known as spectral methods holds much promise. 
For smooth solutions, the errors produced by spectral 
methods decrease exponentially as computational re- 
sources are increased, whereas the errors of finite dif- 
ference methods, the methods used by the majority of 
binary black hole simulations, decrease polynomially. In- 
deed, spectral methods have been used to produce very 
accurate initial data for binary black holes and neutron 
stars [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56], 
and they have been used to produce the longest and 
most accurate binary black hole inspiral simulation to 
date [9, 57]. 

However, a key difficulty with time-dependent spec- 
tral binary black hole simulations has been handling the 
merger of the two holes. For example, the spectral sim- 
ulations described in [9, 12, 57] are very accurate and 
efficient, but they follow only the inspiral of the two 
black holes, and fail just before the holes merge. This 
is sufficient for some applications, such as comparing 
post-Newtonian formulae with numerical results during 
the inspiral and finding accurate analytic templates that 
match the numerical inspiral waveforms [9, 12], but for 
most purposes the merger is the most crucial part of the 
process: for instance the gravitational-wave emission is 
the strongest during merger, and details of the merger 
determine the recoil velocity of the final black hole. 

In this paper we present a spectral binary black hole 
simulation that follows sixteen orbits of the binary plus 
merger and ringdown of the merged black hole. In Sec. II 
we describe the equations, gauge conditions, and numer- 
ical methods we use to solve Einstein's equations; in par- 
ticular, Sees. lie and II D describe changes to our gauge 
conditions that allow simulation of the merger, and our 
method for extending the evolution through ringdown. In 
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Sec. Ill we discuss extraction of the gravitational wave- 
form from the simulation, including the process of extrap- 
olating the waveform to infinity. Sec. Ill also includes an 
estimate of the uncertainty in the waveform from sev- 
eral sources. Finally, in Sec. IV we discuss outstanding 
difficulties and future improvements. 



II. SOLUTION OF EINSTEIN'S EQUATIONS 

A. Initial data 

The initial data describe two nonspinning black holes, 
each with Christodoulou mass M /2, in quasicircular or- 
bit with low eccentricity. The initial data are exactly as 
described in Ref. [9] . Briefly, initial data are constructed 
within the conformal thin sandwich formalism [58, 59] 
using a pscudospectral elliptic solver [49]. We employ 
quasiequilibrium boundary conditions [50, 60] on spher- 
ical excision boundaries, choose conformal flatness and 
maximal slicing, and use Eq. (33a) of Ref. [53] as the 
lapse boundary condition. The spins of the black holes 
are made very small (~ 10^^) via an appropriate choice of 
the tangential shift at the excision surfaces, as described 
in [53] . Finally, the initial orbital eccentricity is tuned to 
a very small value 5 x 10~^) using the iterative proce- 
dure described in Ref. [9], which is an improved version 
of the procedure of Ref. [Gl]. 

B. Evolution of the inspiral phase 

The evolution of the first ~ 15 binary orbits is identical 
to the simulation presented in Ref. [9]. We describe it 
here briefly in order to facilitate the presentation of our 
method for continuing the evolution through merger and 
ringdown, which is described in Sees. II C and II D. 

The Einstein evolution equations are solved with the 
pscudospectral evolution code described in Ref. [57]. 
This code evolves a first-order representation [62] of the 
generalized harmonic system [63, 64, 65]. We handle 
the singularities by excising the black hole interiors from 
the computational domain. Our outer boundary condi- 
tions [62, 66, 67] are designed to prevent the influx of un- 
physical constraint violations [68, 69, 70, 71, 72, 73, 74] 
and undesired incoming gravitational radiation [75, 76], 
while allowing the outgoing gravitational radiation to 
pass freely through the boundary. 

We employ the dual-frame method described in 
Ref. [57]: we solve the equations in an "inertial frame" 
that is asymptotically Minkowski, but our domain de- 
composition is fixed in a "comoving frame" that rotates 
with respect to the inertial frame and also shrinks with 
respect to the inertial frame as the holes approach each 
other. The positions of the holes are fixed in the comov- 
ing frame; we account for the motion of the holes by dy- 
namically adjusting the coordinate mapping between the 
two frames. Note that the comoving frame is referenced 



only internally in the code as a means of treating mov- 
ing holes with a fixed domain. Therefore all coordinate 
quantities (e.g. black hole trajectories, wave-extraction 
radii) mentioned in this paper are inertial-frame values 
unless explicitly stated otherwise. 

As described in [9], the mapping between inertial and 
comoving coordinates for the inspiral, expressed in polar 
coordinates relative to the center of mass of the system, 
is 



,J2 



ait) + (1 - ait)) 



= ^' + bit), 



(1) 

(2) 
(3) 



where a(i) and bit) are functions of time, and Rq is a con- 
stant usually chosen to be roughly the radius of the outer 
boundary in comoving coordinates. Here primes denote 
the comoving coordinates. For the choice R'q = oo, the 
mapping is simply a rotation by bit) plus an overall con- 
traction given by ait). The functions ait) and bit) are 
determined by a dynamical control system as described 
in Ref. [57]. This control system dynamically adjusts 
ait) and bit) so that the centers of the apparent horizons 
remain stationary in the comoving frame. Note that the 
outer boundary of the computational domain is at a fixed 
comoving radius i?maxj so the inertial-coordinate radius 
of the outer boundary -Rinax(0 ^ function of time. 

The gauge freedom in the generalized harmonic system 
is fixed via a freely specifiable gauge source function Ha 
that satisfies the constraint 



^ Ca — Tab^ 



Ha, 



(4) 



where V^bc are the spacetime Christoffel symbols. To 
choose this gauge source function, we first define a new 
quantity Ha that has the following two properties: (1) Ha 
transforms like a tensor, and (2) in inertial coordinates 
Ha = Ha- We choose Ha so that the constraint equa- 
tion (4) is satisfied initially, and we demand that Ha' is 
constant in the moving frame, i.e., that dt'Ha' — 0. 



C. Extending inspiral runs through merger 

If the inspiral runs described above are allowed to con- 
tinue without any modification of the algorithm, then 
as the binary approaches merger, the horizons of the 
black holes become extremely distorted and the dynam- 
ical fields begin to develop sharp (but numerically con- 
vergent) features near each hole. These features grow 
rapidly in time, eventually halting the simulation before 
merger. This is due to a gauge effect: The gauge condi- 
tion used during the inspiral, namely fixing Ha in time 
in the comoving frame, was chosen based on the idea 
that each black hole is in quasiequilibrium in this frame. 
Once the black holes begin to interact strongly, this gauge 
condition no longer allows the coordinates to sufficiently 
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react to the changing geometry, and coordinate singular- 
ities develop. 

Therefore we must modify our gauge conditions in or- 
der to handle merger. Because the inspiral gauge works 
so well before merger, we choose to remain in that gauge 
until some time t = tg, and then we change (smoothly) 
to a new gauge. 

We have experimented with several gauge condi- 
tions [77], but so far the simplest gauge choice that works, 
and the one used in the simulations presented here, is 
based on the gauge treatment of Pretorius [1, 65, 78]: 
We promote the gauge source function Ha to an inde- 
pendent dynamical field that satisfies 



gauge does not change too rapidly. We choose 



W,Ha = Qaix,t,ijab)+^2t''dbHa, 



(5) 



where V^Vc is the curved space scalar wave operator (i.e. 
each component of Ha is evolved as a scalar) , ipab is the 
spacetime metric, and t"" is the timelike unit normal to 
the hypersurface. The driving function Qa is 



Q^ 



9{x, t)^3 



1 - N 
iV2- 



(6) 
(7) 



Here N and A^* are the lapse function and the shift vector, 
77, ^1, ^2, and ^3 are constants, and f{x,t) and g{x,t) 
are prescribed functions of the spacetime coordinates (we 
describe our choices for these objects below). 

Equation (5) is a damped, driven wave equation with 
damping parameter ^2 and driving function Qa- The 
driving term Qt in Eq. (6) was introduced by Preto- 
rius [1, 65] to drive the lapse function towards unity so 
as to prevent it from becoming small. The driving term 
Qi is new; it drives the shift vector towards zero near 
the horizons. This causes the horizons to expand in co- 
ordinate space, and has the effect of smoothing out the 
dynamical fields near the horizon and preventing gauge 
singularities from developing. A different gauge choice 
that causes similar coordinate expansion of the horizons 
was introduced in Rcf. [79]. Care must be taken so that 
the horizons do not expand too quickly relative to the ex- 
cision boundaries; otherwise the characteristic fields will 
fail to be purely outgoing (into the holes) at the exci- 
sion boundaries, and excision will fail. We find that with 
appropriate choices of ^1, ^3, f{x,t), and g{x,t) as de- 
scribed below, the horizons expand gradually and not 
too rapidly. 

For the runs presented here we choose = 4, ^1 = 0.1, 
^2 = 10, and ^3 = 0.4. The functions f{x,t) and g{x,t) 
in Eqs. (6) and (7) are chosen based on two criteria: the 
first is that the driving terms Qa are nonzero only near 
the black holes where they are needed; if these terms are 
nonzero in the wave-extraction zone they lead to compli- 
cated gauge dynamics in this region, making waveform 
extraction difficult. The second criterion is that the driv- 
ing terms are turned on in a gradual manner so that the 



fix,t) = gix.t) - (2 - 



X (1 



/2 / 2 



(8) 



where r' is the coordinate radius in comoving coordi- 
nates, and the constants are ai ~ 17. 5M, tT2 ^ 15A/, 
and 0-3 ~ 40Af. Here M is the sum of the initial 
Christodoulou masses of the two holes. 

Equation (5) is a second-order hyperbolic equation, 
which we evolve in first-order form by defining new fields 
H^ and representing (up to the addition of con- 
straints) the appropriate time and space derivatives of 
Ha, respectively: 

Hf = -t^dbHa, (9) 

'^fa - d,Ha. (10) 

The representation of wave equations of this type in first- 
order form is well understood, see e.g., Refs. [62, 80]; the 
result for Eq. (5) is 



ka ' 



(11) 



dtU^ = N>^dkTl^ -Ng^'Ok'^g-^f^N'^dkHa 

+ ^^N'^^^a + NiV'^j ~ g'''d,N)^t 

+ NKn^ + Qa, (12) 

dt<^fa = N^dk<^fa-Nd,IV^ +^-Nd,Ha 



(13) 



where gij is the spatial metric and K is the trace of the 
extrinsic curvature. We choose the constraint-damping 
parameter 7^^ to be ^2 = 4/M. 

These equations are symmetric hyperbolic, and require 
boundary conditions on all incoming characteristic fields 
at all boundaries. The characteristic fields for Eqs. (11)- 
(13) in the direction of a unit spacelike covector rii are 



jjH± 



Hf ±71^$«_^H^^a, 
Ha, 



(14) 
(15) 
(16) 

The (coordinate) characteristic speeds for Ua^, , 
and Zf^ are ±iV — UiN"^ , 0, and —riiN'^, respectively. 

At the excision boundaries all characteristic fields are 
outgoing (i.e. into the holes) or nonpropagating, so no 
boundary conditions are necessary and none arc imposed. 
At the outer boundary, we must impose boundary con- 
ditions on Ua~ and Z^^'^. Define 

A(C/f ^) ^ 9tnf ± n^dt'Pg - l^dtHa, (17) 
A(^D = dtHa, (18) 
A(^D = (<5f-n,n'=)a,$f^, (19) 

where the time derivatives on the right-hand side are 
evaluated using Eqs. (11)-(13). Then we impose the fol- 
lowing boundary conditions: 



dtZ"' 



H2\ 



(20) 

2nfcA^'=nJ"%$f,. (21) 
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Equation (20) is the outgoing-wave boundary condition 
described in detail in Ref. [80]. Equation (21) ensures 
that violations of the artificial constraint = <i>,^ — 
diHa = do not enter the domain through the bound- 
ary; it is the direct analogue of the constraint-preserving 
boundary condition we apply to the analogous variable in 
the generalized harmonic formulation of Einstein's equa- 
tions, Eq. (65) of Ref. [G2]. 

Note that Eqs. (11)-(13) involve only first derivatives 
of the spacetime metric, and similarly, the generalized 
harmonic Einstein equations involve only first derivatives 
of Ha- Therefore, adding Eqs. (11)-(13) to the system 
does not change the hyperbolicity or characteristic fields 
of the generalized harmonic Einstein equations, so we can 
impose the same boundary conditions on the generalized 
harmonic variables as we do during the inspiral, as de- 
scribed in Rcfs. [57, 66]. 

Equations (11)-(13) require as initial data the values of 
Ha and 11^ at t = tg. These quantities can be computed 
from the gauge choice used during the inspiral for t < tg, 
so we choose them to be continuous at t = tg. 

Note that Eqs. (11)-(13) and the boundary condi- 
tions (20) and (21) are written in the inertial coordinate 
system. The equations are actually solved in the comov- 
ing coordinate system using the dual-frame method de- 
scribed in Ref. [57]. 

With the modifications to the gauge conditions de- 
scribed here, the evolution of the binary can be tracked 
up until (and shortly after) the formation of a common 
horizon that encompasses both black holes. Because of 
the more rapid dynamics and the distortions of the hori- 
zons during the merger, we typically increase the numer- 
ical resolution slightly when we make these changes to 
the gauge conditions (this is the difference between the 
first and second entry in the iVpts column in Table I). 
After the common horizon forms, the problem reduces to 
evolving a single highly distorted dynamical black hole, 
rather than two separate black holes. We change the al- 
gorithm to take advantage of this, as described in the 
next section. 



D. Evolution from merger through ringdown 

We make three main changes to our evolution algo- 
rithm once we detect a common apparent horizon. First, 
because there is now only one black hole and not two, 
we interpolate all variables onto a new computational 
domain that contains only a single excised region. Sec- 
ond, we choose a new comoving coordinate system (and 
a corresponding mapping to inertial coordinates) so that 
the new excision boundary tracks the shape of the (dis- 
torted, rotating, pulsating) apparent horizon in the in- 
ertial frame, and so that the outer boundary behaves 
smoothly in time. Third, we modify the gauge condi- 
tions so that the shift vector is no longer driven towards 
zero, allowing the solution to eventually relax to a time- 
independent state. We now discuss these three changes 



in detail. 

Our new computational domain contains only a single 
excised region, and is much simpler than the one used 
until merger. It consists only of nested spherical-shell 
subdomains that extend from a new excision boundary 
^min' chosen to be slightly inside the common apparent 
horizon, to an outer boundary -R^ax that coincides with 
the outer boundary of the old domain. 

To understand how we choose our new comoving frame, 
first recall that in the dual- frame technique [57] , the co- 
moving frame is the one in which the computational do- 
main is fixed, the inertial frame is the one in which the 
coordinates are Minkowski-like at infinity, and the two 
frames are related by a mapping that is chosen so that 
the computational domain tracks the motion of the black 
holes. Let a;° represent the inertial coordinates (which 
are the same before and after merger), let x"^ represent 
the old comoving coordinates, and let x"" represent the 
new comoving coordinates. The mapping between a;'" 
and is given by Eqs. (l)-(3). The mapping between 
x"" and x°- is chosen to be 



l + sin2(7rf/2<aJ 



R" 



+ il-Ait)) 



R'' 



^•'inax-^'-O 



(22) 



t=0 m=~l 



= r"-q{r")Y, ^ \i^{t)Y,„,{e" A"), (23) 
Bit), 



9", 

J," 



(24) 
(25) 



where R'^ax the outer boundary of the premerger com- 
putational domain in the old comoving coordinates, and 
q{r"), A{t), B{t), and Xemit) are functions we will now 
discuss. 

First we describe the angular map: The function B{t) 
is chosen so that the new comoving frame initially rotates 
with respect to the inertial frame, but this rotation slows 
to a halt after a short time. In particular. 



B{t) =Bo + {Bi + B2{t - t,„))e-(*-*")/^«. 



(26) 



where the constants Bq, Bi, and B2 are chosen so that 
B{t) matches smoothly onto b{t) from Eq. (3): B{tm) ~ 
b{tm), B{t,n) = h{tm), and B{t.„,) = b{tm)- Here t^ 
is the time at which we transition to the new domain 
decomposition. The constant Tg is chosen to be on the 
order of 20 M. 

The radial map is a composition of two individual 
maps: Eqs. (22) and (23). The purpose of Eq. (22) 
is to match the outer boundary of the new domain 
smoothly onto that of the old domain, while far from 
the outer boundary Eq. (22) approaches the identity. 
We have found that without the use of Eq. (22), the 
(inertial-coordinate) location of the boundary changes 
nonsmoothly at t = , thereby generating a spurious in- 
going gauge pulse that spoils waveform extraction. The 
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function A{t) is 

A{t) =Ao + (A, + A^it. - ^„0)e-^*-*'"^/"^ (27) 

where the constants Ag, Ai, and A2 are chosen so that 
A{t) matches smoothly onto a{t) from Eq. (1): A{tm) = 
a(t„i), A{tjn) = a(tm), and A{tm) = a{tm)- The constant 
TA is chosen to be on the order of 5M. 

The other piece of the radial map, Eq. (23), is chosen 
so that the apparent horizon is nearly spherical in the 
new comoving coordinates a;"°. The fmiction q{r") is 

g(r") = e-('-"-<H)V--\ (28) 

where i?^H radius of the apparent horizon in co- 

moving coordinates, and aq is a constant of order 20M. 
This function q{r") ensures that the piece of the radial 
map represented by Eq. (23) acts only in the vicinity of 
the merged hole and not in the exterior wave-extraction 
region. 

We now discuss the choice of the functions Xim (t) that 
appear in Eq. (23). Given the known location of the 
apparent horizon in inertial coordinates, the Xim{t) de- 
termine the shape of the apparent horizon in comoving 
coordinates. At t = i^, we choose these quantities so 
that the apparent horizon is spherical (up to spherical 
harmonic component i = f max) in comoving coordinates: 
that is, if the comoving-coordinate radius of the apparent 
horizon as a function of angles is written as 



''AH I 



E 

^=0 m= 



E 



(29) 



then for 1 < < £max we choose Xem{tm) so that 
Qim{tm) = 0. In addition, we choose Aoo(im) = 0; this 
determines i^Xn- ^ > Xgmit) are determined 

by a dynamical feedback control system identical to the 
one described in Ref. [57], which adjusts these functions 
so that the apparent horizon is driven to a sphere (up to 
spherical harmonic component £ = ^max) in comoving co- 
ordinates. This dynamical feedback control allows us to 
freely choose the first and second time derivatives of Xim 
ai t = tm- Simply choosing these to be zero causes the 
control system to oscillate wildly before settling down, 
and unless the time step is very small, these oscillations 
are large enough that the excision boundary crosses the 
horizon and our excision algorithm fails. So instead, we 
obtain the time derivatives of Xim by finding the apparent 
horizon at several times surrounding t = tm, comput- 
ing Afin at these times, and finite-differencing in time. 
For the equal-mass zero-spin merger presented here, in 
Eq. (23) it suffices to sum only over even £ and m and to 

choose £max = 6. 

The last change we make before continuing the simu- 
lation past merger is to modify the functions f{x, t) and 
g{x,t), which before merger were given by Eq. (8), to 



(2 

X (1 ~ e-^*-*«^'/'"2')e- 
ff(x,t) = /(x,t)e-(*-*")^/'^S 



-(t-tg)/0-l^ 



2 



Run 


R' 

-''-max 


R" 


R'o Npts 


CPU-h CPU-h/T 


30cl/N4 


462 


462 


698 (57^59^57^) 


8,800 


2.0 


30cl/N5 


462 


462 


698 (62^66^63^) 


15,000 


3.4 


30cl/N6 


462 


462 


698 (67^73^70^) 


23,000 


5.3 


30c2/N6 


722 


96 


oo (71^76^63^) 


25,000 


5.7 



(30) 
(31) 



TABLE I: Outer boundary parameters, collocation points, 
and CPU usage for several zero-spin binary black hole evo- 
lutions. The first column identifies the inspiral run in the 
nomenclature of Ref. [9]. Apts is the approximate number 
of collocation points used to cover the entire computational 
domain. The three values for Apts are those for the inspiral, 
merger, and ringdown portions of the simulation, which are 
described in Sections II B, II C, and II D, respectively. The 
outer boundary parameters -Rmax, -Rmax and R'o, as well as 
run times T, are in units of the initial Christodoulou mass M 
of the system, which provides a natural time and length scale. 



where (74 ~ 7M. The modification of g{x,t) turns off 
the term in the gauge evolution equations that drives 
the shift to zero near the holes. Before merger, it is ad- 
vantageous to have the shift driven to zero so that the 
horizons expand in coordinate space and so that growing 
gauge modes remain inside the common horizon. After 
merger, however, it is no longer desirable for the horizon 
to expand, since this would prevent the solution from 
eventually settling down to a time-independent state in 
which the horizon is stationary with respect to the coor- 
dinates. 

To summarize, the steps involved in the transition from 
evolving a binary black hole spacetime to evolving a 
merged single black hole spacetime are as follows: (1) 
Find the common apparent horizon in the inertial frame 
at several times near t — t^- (2) Solve for the A^,„(im) 
that make the horizon spherical in the comoving frame, 
and simultaneously solve for i?AH- (3) Choose the in- 
ner boundary of the new computational domain R'^-^^^ to 
be slightly less than i^^H' ^'^'^ choose the outer bound- 
ary i?"jax [for sufficiently small a(tm) it is necessary to 
choose i?max < R'max that the mapping (22) is invert- 
ible]. At this point the computational domain and the 
mapping (22)- (25) have been determined. (4) Interpo- 
late all dynamical variables from the old computational 
domain onto the new one. This interpolation is done via 
the spectral expansion in the old domain, so it introduces 
no additional error. (5) Modify the gauge source evolu- 
tion equations so that the shift is no longer driven to 
zero. (6) Continue the evolution on the new computa- 
tional domain. All of these steps can be automated. 



E. Properties of the numerical solution 

In Table I we list outer boundary parameters, resolu- 
tions, and run times of several runs we have done using 
the algorithm described above. Three of these runs are 
identical except for numerical resolution, and the fourth 
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FIG. 1: Spacetime diagram showing the spacetime volume 
simulated by the numerical evolutions listed in Table I. Each 
curve represents the worldline of the outer boundary for a 
particular simulation. The magnified views on the right show 
that the outer boundary moves smoothly near merger. The 



transition times tg = 
on the right panels. 



3917Af and t„ 



3940M are indicated 



FIG. 2: Constraint violations of run 30c 1. The top panel 
shows the norm of all constraints, normalized by the 
norm of the spatial gradients of all dynamical fields. The 
bottom panel shows the same data, but without the normal- 
ization factor. The norms are taken over the portion of the 
computational volume that lies outside apparent horizons. 



is performed on a different domain with a different outer 
boundary location. As discussed above, the outer bound- 
ary of our simulation varies in time because of the dual- 
frame approach we use to follow the black holes. Figure 1 
is a spacetime diagram illustrating the region of space- 
time being evolved in our simulation. 

We do not explicitly enforce either the Einstein con- 
straints or the secondary constraints that arise from writ- 
ing the system in first-order form. Therefore, examin- 
ing how well these constraints are satisfied provides a 
useful consistency check. Figure 2 shows the constraint 
violations for run 30cl. The top panel shows the 
norm of all the constraint fields of our first-order gener- 
alized harmonic system, normalized by the norm of 
the spatial gradients of the dynamical fields (see Eq. (71) 
of Ref. [G2]). The bottom panel shows the same quan- 
tity, but without the normalization factor [i.e., just the 
numerator of Eq. (71) of Ref. [62]]. The norms are 
taken over the portion of the computational volume that 
lies outside apparent horizons. At early times, t < 500M, 
the constraints converge rather slowly with resolution be- 
cause the junk radiation contains high frequencies. Con- 
vergence is more rapid during the smooth inspiral phase, 
after the junk radiation has exited through the outer 
boundary. 

The constraints increase as the holes approach each 
other and the solution becomes increasingly distorted. 



At t = 3917M {t = 3927Af for resolution N4), the gauge 
conditions are changed (cf. Sec. II C) and the resolution 
is increased slightly (compare the first and second entry 
in the A^pts column in Table I). Because of the change 
of resolution, the constraints drop rapidly by almost two 
orders of magnitude, but then they begin to grow again. 
The transition to a single-hole evolution (cf. Sec. II D) 
occurs at t = 3940M {t = 3948M for resolution N4). 
At this time the constraint norm drops by about an or- 
der of magnitude because the region in which the largest 
constraint violations occur — the interior of the common 
horizon — is newly excised. 

After the binary proceeds through inspiral, merger, 
and ringdown, it settles down to a final stationary black 
hole. In our simulation this final state is not expressed 
in any standard coordinate system used to describe Kerr 
spacetime, but nevertheless the final mass and spin of 
the hole can be determined. The area A of the apparent 
horizon provides the irreducible mass of the final black 
hole, 

Mi„ = y^A/16n, (32) 

which we find to be Mi„/M = 0.88433 ± 0.00001, where 
M is the sum of the initial irreducible masses of the black 
holes. The uncertainty in Mi„/M is determined from the 
difference between runs 30cl/N6, 30cl/N5, and 30c2/N6, 
so it includes only uncertainties due to numerical reso- 
lution and outer boundary location. We have verified 



7 



Initial orbital eccentricity: 


e 


~ 5 X 10"^ 


Initial spin of each hole: 






Time of evolution: 


T/M 


= 4330 


Final Christodoulou mass: 


Mf/M 


= 0.95162 ± 0.00002 


Final spin: 




= 0.68646 ± 0.00004 



TABLE II: Physical parameters describing the equal-mass 
nonspinning binary black hole evolutions presented here. 
The dimensionful quantity M is the initial sum of the 
Christodoulou masses of the black holes. Uncertainty esti- 
mates include numerical uncertainties and the effects of vary- 
ing the outer boundary location. 



that the uncertainty due to the finite resolution of our 
apparent horizon finder is neghgible. 

The final spin Sf of the black hole can be computed 
by integrating a quasilocal angular momentum density 
over the final apparent horizon [81, 82]. Our imple- 
mentation of this method is described in detail in Ap- 
pendix A of [5-5]. Furthermore, an alternative method of 
computing the final spin, which is based on evaluating 
the extremal values of the 2-dimensional scalar curva- 
ture on the apparent horizon and comparing these values 
to those obtained analytically for a Kerr black hole, is 
also described in [55]. Using these measures, we deter- 
mine the dimensionless spin of the final black hole to 
be Sf/Mj = 0.68646 ± 0.00004, where the uncertainty 
is dominated by the difference between runs 30cl/N6 
and 30cl/N5 rather than by the differences between dif- 
ferent methods of measuring the spin. Here Mf is the 
Christodoulou mass of the final black hole. 



M 



f 



5^2 



ami; 



(33) 



We find that the ratio of the final to initial black hole 
mass is Mf/M = 0.95162 ± 0.00002. The mass and spin 
of the final hole are consistent with those found by other 
groups [2, 4, 83, 84, 85]. Physical parameters describing 
the evolutions are summarized in Table II. 



III. COMPUTATION OF THE WAVEFORM 



Penrose scalar ^^4, following the same procedure as in 
Refs. [61, 86]. To summarize, we compute 



*4 = -Cc.^,pJ^'tm''mP 



where 



V2 



1/9 I dY 
i- 



V2r \d6 sin 6190 



(34) 

(35a) 
(35b) 



Here (r, 9, </>) denote the standard spherical coordinates 
in the inertial frame, is the timelike unit normal to the 
spatial hypersurface, and is the outward-pointing unit 
normal to the extraction sphere. We then expand ^'4 in 
terms of spin- weighted spherical harmonics of weight —2: 



Y,^'r{t.r)-2YUe,c^), (36) 



where the 'i'4" are expansion coefficients defined by this 
equation. 

Note that our choice of is not exactly null nor 
exactly of unit magnitude at finite r, as is required by 
the standard definition. The resulting ^'4'" computed at 
finite r will therefore disagree with the waveforms ob- 
served at infinity. Our definition does, however, agree 
with the standard definition of ^'4™ as r — > 00. Because 
we extrapolate the extracted waves to find the asymp- 
totic radiation field (see Sec. IIIC), these tetrad effects 
should not play a role: Relative errors in ^'4™ introduced 
by using the simple coordinate tetrad fall off like pow- 
ers of M/r, and thus should vanish after extrapolating 
to obtain the asymptotic behavior. More careful treat- 
ment of the extraction method — such as those discussed 
in Refs. [87, 88, 89] — may improve the quality of extrap- 
olation and would be interesting to explore in the future. 

In this paper, we focus on the dominant (i, m) = (2, 2) 
mode. Following common practice (see e.g. [84, 85]), we 
split the extracted waveform into real phase 4> and real 
amplitude A, defined by 



The numerical solution of Einstein's equations ob- 
tained using the methods described above yields the 
spacetime metric and its first derivatives at all points 
in the computational domain. In this section we describe 
how this solution is used to compute the key quantity 
relevant for gravitational-wave observations: the gravi- 
tational waveform as seen by an observer infinitely far 
from the source. 



A. Waveform extraction 

Gravitational waves are extracted from the simulation 
on a sphere of coordinate radius r using the Newman- 



ir,t) 



.4^(r,t) =^(r,Oe- 
The gravitational-wave frequency is given by 



dt 



(37) 



(38) 



The minus sign in Eq. (37) is chosen so that the phase 
increases in time and to is positive. 

The {I, m) ~ (2, 2) waveform, extracted at a single 
radius for nm 30cl/N6, is shown in Fig. 3. The short 
pulse at t ~ 200A/ is caused by imperfect initial data 
that are not precisely in equilibrium; this pulse is usually 
referred to as "junk radiation" . 
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FIG. 3: Gravitational waveform extracted at finite radius r — 
225M, for the case 30cl/N6 in Table I. The left panel zooms 
in on the inspiral waveform, and the right panel zooms in on 
the merger and ringdown. 



B. Convergence of extracted waveforms 

In this section we examine the convergence of the grav- 
itational waveforms extracted at fixed radius, without 
extrapolation to infinity. This allows us to study the 
behavior of our code without the complications of ex- 
trapolation. The extrapolation process and the resulting 
extrapolated waveforms are discussed in Sec. IIIC. 

Figure 4 shows the convergence of the gravitational- 
wave phase (f> and amplitude A with numerical resolu- 
tion. For this plot, the waveform was extracted at a fixed 
inertial-coordinate radius of r = 60M. This fairly small 
extraction radius was chosen to allow a comparison of 
the simulations 30cl and 30c2. Each solid line in the top 
panel shows the absolute difference between </> computed 
at some particular resolution and computed from our 
highest-resolution run, labeled 30cl/N6 in Table I. The 
solid curves in the bottom panel similarly show the rela- 
tive amplitude differences. When subtracting results at 
different resolutions, no time or phase adjustment has 
been performed. The noise at early times is due to "junk 
radiation" generated near t = 0. While most of this radi- 
ation leaves through the outer boundary after one cross- 
ing time, some remains visible for a few crossing times^. 
The plots show that the phase difference accumulated 
over 16 orbits plus merger and ringdown is less than 0.1 
radians for our medium resolution, and the relative am- 
plitude differences are less than 0.015; these numbers can 
be taken as an estimate of the numerical truncation error 
of our medium resolution run. 

Also shown as a dotted curve in each panel of Fig. 4 
is the difference between our highest-resolution run, 
30cl/N6, and a similar run but with a different outer 
boundary location, 30c2/N6. The 30c2 run initially has 
a more distant outer boundary than 30cl, but during the 
inspiral the outer boundary moves rapidly inward, as seen 
in Fig. 1, so that extraction of the full waveform is pos- 



^ The junk radiation at early times is discussed in more detail 
in Ref. [9] (specifically, just before Eq. (9) and in the third 
paragraph of Sec. II E), which presents the exact same waveform 
as shown here but without merger and ringdown. 





1000 2000 3000 4000 
t/M 



FIG. 4: Convergence of waveforms with numerical resolution 
and outer boundary location. Shown are phase and amplitude 
difi'erences between numerical waveforms ^^l^ computed using 
different numerical resolutions. Shown also is the difference 
between our highest-resolution waveforms using two different 
outer boundary locations. All waveforms are extracted at 
r = 60M, and no time shifting or phase shifting is done to 
align waveforms. 



sible only for extraction radii r < 75 M . Comparing runs 
30cl and 30c2 provides an estimate of the uncertainty 
in the waveform due to outer boundary effects such as 
imperfect boundary conditions that might reflect outgo- 
ing waves. From Fig. 4 we estimate this uncertainty to 
be 0.03 radians in phase and half a percent in amplitude 
(when no time shift is applied). 

Figure 5 is the same as Fig. 4 except each waveform is 
time shifted and phase shifted so that the maximum am- 
plitude of the wave occurs at the same time and phase. 
This type of comparison is relevant for analysis of data 
from gravitational-wave detectors: when comparing ex- 
perimental data with numerical detection templates, the 
template will be shifted in both time and phase to best 
match the data. For this type of comparison. Fig. 5 shows 
that the numerical truncation error of our medium resolu- 
tion run is less than 0.01 radians in phase and 0.1 percent 
in amplitude for t > lOOOAf. At earlier times, the errors 
are somewhat larger and are dominated by residual junk 
radiation. Our uncertainty due to outer boundary effects 
is similar to that in Fig. 4: about 0.02 radians in phase 
and half a percent in amplitude. Boundary effects are 
most prominent during the ringdown. 
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and the tortoise-coordinate radius [90] is 



2000 3000 4000 
t/M 



FIG. 5: Convergence of waveforms with numerical resolution 
and outer boundary location. Same as Fig. 4 except wave- 
forms are time-shifted and phase-shifted so that the maximum 
amplitude occurs at the same time and phase. 



C. Extrapolation of waveforms to infinity 

Our numerical simulations cover only a finite spacetime 
volume, as shown in Fig. 1, so it is necessary to extract 
our numerical waveforms at a finite distance from the 
source. However, gravitational-wave detectors measure 
waveforms as seen by an observer infinitely far from the 
source. Accordingly, after extracting waveforms at multi- 
ple finite radii, we extrapolate these waveforms to infinite 
radius using a procedure similar to that described in [9]. 
This extrapolation procedure is intended to remove not 
only near-field effects that are absent at infinity, but also 
gauge effects that can be caused by the time-dependence 
of the lapse function or the nonoptimal choice of tetrad 
for computing ^'4. 

The extraction procedure described in Sec. Ill A yields 
a set of waveforms ^|^(t,r), with each waveform ex- 
tracted at a different radius. To extrapolate to infinite 
radius we must compare waveforms at different radii, but 
these waveforms must be offset in time by the light-travel 
time between adjacent radii. To account for this time 
shift, for each extraction radius we compute ^'4^(it,r), 
where u is the retarded time at that radius. Assuming 
for simplicity that the background spacetime is nearly 
Schwarzschild, we compute the retarded time u using 

u = t,-r*, (39) 

where t,, is some approximation of Schwarzschild time. 



^arcal + 2A/aDM In 



^arcal 



(2A/A 



DM 



(40) 



Here Madm is the ADM mass of the initial data, 
and Taroai = \/ A/in, where A is the measured (time- 
dependent) area of the extraction sphere. If we were 
to choose ts to be simply the coordinate time t, then 
the retarded time coordinate u would fail to be null, 
largely because the lapse function in our simulation is 
time-dependent and differs from the Schwarzschild value. 
We attempt to account for this by assuming that our 
background spacetime coordinates are Schwarzschild, but 
with gtt replaced by — iVj^,g, where iVavg is the (time- 
dependent) average value of the lapse function measured 
on the extraction sphere. Under these assumptions, it 
can be shown that the one-form 



avg 



2MADM/?'arcal 



dt - dr* 



(41) 



is null, so we equate this one-form with du and thus define 



ts = 



^1 — 2MADM/''aroal 



-.dt. 



(42) 



We show below (cf. Fig. 9) that choosing Eq. (42) in- 
stead of is = t significantly increases the accuracy of our 
extrapolation procedure during merger and ringdown. 

Having computed the retarded time at each extrac- 
tion radius, we now consider the extracted waveforms as 
functions of retarded time u and extraction radius raroah 
i.e. '^4^{u, Tareai). At each value of u, we have the phase 
and amplitude of ^"1^ at several extraction radii raroai- 
Therefore at each value of u, we fit phase and amplitude 
separately to a polynomial in 1 /raroai: 



raroai) = <^(0) (") + X! T^' ^ ^ 
fe=l 
n 

raroai A{u, r) = A(o) (u) + ^ 



1^—1 ^arcal 

Aik){u) 



(43) 
(44) 



fc=i 



roal 



The phase and amplitude of the desired asymptotic wave- 
form are thus given by the leading-order term of the ap- 
propriate polynomial, as a function of retarded time: 



= '/'(o)(u), 

raroai ^(m) = ^(0)(w)- 



(45) 
(46) 



Figure 6 shows phase and amplitude differences be- 
tween extrapolated waveforms that are computed us- 
ing different values of polynomial order n in Eqs. (43) 
and (44). For the extrapolation we use waveforms ex- 
tracted at radn 75M, 85M, lOOM, llOM, 130M, 140M, 
150A'f, 160Af, 170M, 180M, 190M, 200M, 210M, and 
225M. From Fig. 6 it is clear that increasing n increases 
the accuracy of the extrapolation in smooth regions, but 
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FIG. 6: Convergence of extrapolation to infinity for extrap- 
olation of order n. For each n, plotted is the extrapolated 
waveform from run 30cl/N6 using order n + 1 minus the ex- 
trapolated waveform using order n. The top panel shows 
phase difi'erences, the bottom panel shows amplitude difi'er- 
ences. No shifting in time or phase has been done for this 
comparison. Increasing n increases accuracy in smooth re- 
gions but also amplifies noise. 
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FIG. 7: Late-time phase convergence of extrapolation to in- 
finity. Same as the top panel of Fig. 6, except zoomed to 
late times. The peak amplitude of the waveform occurs at 
ts-r* ^ 3954M. 
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also amplifies any noise present in the waveform. Our 
preferred choice, n — 3, gives a phase error of 0.005 radi- 
ans and a relative amplitude error of 0.003 during most 
of the inspiral, and a phase error of 0.01 radians and a 
relative amplitude error of 0.01 in the ringdown. The 
junk radiation epoch ts — r* % lOOOM has moderately 
larger errors than the ringdown. If we were to choose 
instead n 4, we would gain higher accuracy in the 
smooth regions at the expense of increased noise in the 
junk radiation epoch and slightly larger errors during the 
merger and ringdown. 

Figure 7 is the same as the top panel of Fig. 6, ex- 
cept zoomed to late times. Note that during merger 
and ringdown, the extrapolation procedure does not con- 
verge with increasing extrapolation order n: the phase 
differences are slightly larger for larger n. This lack of 
convergence suggests that the nonextrapolated numerical 
waveform contains some small contamination that does 
not obey the fitting formulae, Eqs. (43) and (44). Fig- 
ure 8 shows the n=:l and n=2 convergence curves from 
Fig. 7, but computed for two different numerical resolu- 
tions, 30cl/N5 and 30cl/N6. The N5 and N6 Hues are 
very close to each other in this figure, indicating that 
the lack of convergence with extrapolation order n is not 
dominated by insufhcient numerical resolution. We sus- 
pect that the main contribution is instead due to gauge 



FIG. 8: Effect of numerical resolution on extrapolation to in- 
finity. The solid curves are identical to the "n=l" and "n=2" 
curves from Fig. 7. The dotted curves are the same quantities 
computed using the lower resolution run 30cl/N5. 



effects. Such gauge effects might be reduced by improv- 
ing the gauge conditions in the numerical simulation or 
by adopting more sophisticated wave extraction and ex- 
trapolation algorithms that better compensate for dy- 
namically varying gauge fields. 

Indeed, we have already made a first attempt at cor- 
recting for a time-dependent lapse function by using tg 
from Eq. (42) to compute the retarded time. Figure 9 
illustrates the importance of this correction. Figures 7 
and 9 differ only in the choice of used to compute the 
retarded time: In Fig. 7, tg is obtained from Eq. (42), 
and in Fig. 9, ts is simply the coordinate time t. Us- 
ing the naive choice ts = t clearly results in much larger 
phase differences that diverge with increasing n and grow 
in time. 

In Fig. 10 we examine the difference between extrapo- 
lated waveforms and waveforms that have been extracted 
at a finite radius. We compare our preferred wave- 
form, 30cl/N6 extrapolated to infinity using n = 3, ver- 
sus nonextrapolated waveforms and versus extrapolated 
waveforms with different values of n. Because the extrap- 
olated and nonextrapolated waveforms differ by overall 
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FIG. 9: Effect of ts on extrapolation to infinity. Same as 
Fig. 7, except the quantity ts that appears in the retarded 
time, Eq. (39), is chosen to be coordinate time t rather than 
the integral in Eq. (42). Note the difference in vertical scale 
between this figure and Fig. 7. 
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FIG. 10: Comparison of extrapolated and nonextrapolated 
waveforms. Plotted are differences between selected wave- 
forms and the 30cl/N6 waveform extrapolated to infinity us- 
ing n = 3. Each selected waveform is labeled by the nu- 
merical resolution (N4, N5, or N6), and either the extraction 
radius (for nonextrapolated waveforms) or the extrapolation 
order (for extrapolated waveforms). Each waveform has been 
shifted in time and phase so as to minimize the least-squares 
difference from the N6, n = 3 waveform. The top panel shows 
phase differences, the bottom panel shows amplitude differ- 
ences. Differences between extrapolated and nonextrapolated 
waveforms are much larger than differences between different 
extrapolation orders. Phase differences between resolutions 
N5 and N6, and amplitude differences between all three reso- 
lutions, are indistinguishable on the plot. 



time and phase offsets which are irrelevant for many pur- 
poses, each waveform in Fig. 10 has been shifted in time 
and phase so as to best match with the n ~ 3 extrap- 
olated waveform. This best match is determined by a 
simple least-squares procedure: we minimize the func- 
tion 



/(io,0o) = 5](Ai(t,) 



J4>iiti 



-A2it,+to)e 



^{4>■2{t^+to)+<t>o) 



(47) 

by varying to and (f)Q. Here Ai, (pi, A2, and 02 are the am- 
plitudes and phases of the two waveforms being matched, 
and the sum goes over all times ti at which waveform 1 
is sampled. 

We find from Fig. 10 that extrapolation to infinity has 
a large effect on the phase of the final waveform and 
a much smaller effect on the amplitude, when compar- 
ing to data extracted at our outermost extraction radius, 
r = 225M. The r = 225M waveforms have an accu- 
mulated phase difference of 0.2 radians relative to the 
extrapolated waveform, much larger than the difference 
between different extrapolation orders or different numer- 
ical resolutions. For extraction at smaller radii, the dif- 
ferences are larger still, the r = 60M waveform having a 
phase difference of 0.8 radians and amplitude difference of 
20 percent compared to the extrapolated waveform. We 
find that the phase differences between extrapolated and 
nonextrapolated waveforms scale quite accuratly like 1 /r, 
and the amplitude differences scale roughly like l/r^-^, 
where r is the extraction radius. These scalings seem 
to be related to near-field effects, for which one expects 
scalings like 1/r in phase and 1/r^ in amplitude [86]. 

Figure 11 presents the final waveform after extrapola- 
tion to infinite radius. There are 33 gravitational- wave 
cycles before the maximum of |*I'4|. The simulation is 
further able to resolve 10 gravitational- wave cycles dur- 
ing ringdown, during which the amplitude \'^4\ drops by 
four orders of magnitude. 



IV. DISCUSSION 

We have presented the first spectral computation of 
a binary black hole inspiral, merger, and ringdown, and 
we have extracted accurate gravitational waveforms from 
our simulation. A key ingredient in handling the merger 
phase is a choice of gauge that causes the individual 
holes to expand in coordinate size. This eliminates the 
coordinate singularities that prevented our earlier sim- 
ulations from continuing through merger. The largest 
downside to the gauge used here is that the success of 
the method depends sensitively on some of the gauge pa- 
rameters, namely cti and (T2 m Eq. (8), and ^1 and ^3 in 
Eqs. (6) and (7). If these parameters are chosen poorly, 
the characteristic fields at the excision boundaries fail to 
be purely outgoing (i.e. into the holes) at some instant 
in time, causing the code to terminate due to lack of a 
proper boundary condition at an excision boundary. An 



12 




1000 



2000 



(tg-r*)/M 



3000 



3800 



3900 



4000 



(t^-r*)/M 



0.003- 
0.001 r 
0.0003 
0.0001 







ReCrM'P^) 1=2, m=2 
I r M ¥4 I 1=2, m=2 



1000 



2000 



3000 



(t^-r*)/M 



3800 



3900 



4000 



(tg-r*)/M 



FIG. 11: Final waveform, extrapolated to infinity. The top panels show the real part of ^4 with a linear y-axis, the bottom 
panels with a logarithmic y-axis. The right panels show an enlargement of merger and ringdown. 



alternative approach to gauge conditions for the general- 
ized harmonic system [77] is in progress, and promises to 
be more robust. 

We compute the spin of the final black hole with three 
distinct diagnostics, one based on approximate rotational 
Killing vectors, the others based on the minimum and 
maximum of the scalar curvature of the apparent horizon 
(XAKV, X™i and X^c^ ^^'^ language of Appendices A 
and B of [55]). We find that these diagnostics agree to an 
exquisite degree. Since these diagnostics coincide exactly 
for a Kerr black hole, this suggests that the final state 
is indeed a Kerr black hole. The uncertainty of the final 
spin quoted in Sec. HE is due to numerical truncation 
error, (i.e. differences between resolutions 30cl/N5 and 
30cl/N6), rather than due to differences between spin 
diagnostics, and we find Sf/Mj = 0.68646 ± 0.00004, 
and Mf = (0.95162 ± 0.00002)M. 

The physical waveform at infinity produced by any nu- 
merical relativity code should of course be independent 
of the coordinates used during the simulation. However, 
in practice it is difficult to remove coordinate effects from 
the waveform for several reasons. First, waveforms are 
typically extracted on coordinate spheres (not geomet- 
ric spheres) of finite radius as functions of coordinate 
time (which may not agree with proper time at infin- 
ity). Second, the extracted waveform on a given sphere 
is typically expanded in spin-weighted spherical harmon- 



ics sY£m{d,(f>) using the 9 and (f> coordinates from the 
simulation rather than some geometrically defined and 
(j) coordinates. Finally, standard formulae equating ^'4 
with the asymptotic radiation field assume that \l/4 is 
computed at infinity. Such gauge ambiguities can be sig- 
nificant for the accuracy of waveforms from numerical 
simulations [87, 88, 89]. Indeed, if we choose a deliber- 
ately "bad" gauge just after merger by omitting the fac- 
tor 1"^ in the function f{x,t) [cf. Equation (30)], 
we find that the lapse function oscillates in time even 
at large distances, and that the resulting waveform ex- 
tracted at a finite radius differs by more than a radian 
in phase from the waveform presented here. We defer 
further discussion of gauge effects on the waveform to a 
future paper. 

We have also shown that extrapolation of waveforms to 
infinity is crucial: waveforms extracted at a finite radius 
differ (particularly in phase) from waveforms extrapo- 
lated to infinity by far more than the numerical errors, 
as shown in Fig. 10. Although it is likely that the need for 
extrapolation may be somewhat reduced by more sophis- 
ticated algorithms for wave extraction at finite radius, it 
appears that most of the difference between waveforms 
that have and have not been extrapolated to infinity is 
due to physics (in the form of near-zone effects) rather 
than to gauge and tetrad ambiguities [86]. 

We are currently extending our methods to binary 
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black holes with unequal masses and nontrivial spins. 
Inspiral simulations for these more generic systems have 
already been computed by our code; it remains to be 
seen whether mergers of more generic black hole systems 
can be simulated with the methods described here, or 
whether alternative gauge conditions, such as those de- 
scribed in Ref. [77], will be necessary. 

It would be interesting to compare the waveforms pre- 
sented here with those from other groups computing bi- 
nary black hole mergers, particularly since other groups 
use different numerical methods, different formulations 
of the equations, and different gauge conditions than our 
group. Several such comparisons are presently under way. 

Waveforms are available at 
http: / /www. black-holes.org/Waveforms. html. 
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